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Abstract 

We compute the complete 0(a 2 ) QED corrections to the electron energy spectrum in unpolarized 
muon decay, including the full dependence on the electron mass. Our calculation reduces the 
theoretical uncertainty on the electron energy spectrum well below 10 -4 , the precision anticipated 
by the TWIST experiment at TRIUMF, which is currently performing this measurement. For 
this calculation, we extend techniques we have recently developed for performing next-to-next-to- 
leading order computations to handle the decay spectra of massive particles. Such an extension 
enables further applications to precision predictions for b, t, and Higgs differential decay rates. 
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I. INTRODUCTION 



The decay of a muon into an electron and a pair of neutrinos, fi — > ez/ M z/ e , occupies an 
important role in particle physics. The measurement of the muon lifetime [jj leads to the 
most accurate determination of the Fermi coupling constant, Gp. The muon anomalous 
magnetic moment is one of the most precisely measured quantities in nature 0, HI, and 
provides important constraints on physics beyond the Standard Model (SM) 0]. Searches 
for lepton flavor-violating decays of the muon, such as /i — > ej and /i — > eee, constrain the 
flavor sector of many SM extensions 

The calculations of radiative corrections to muon decay have a long and storied history 
The one-loop QED corrections were first performed within the Fermi theory of weak interac- 
tions in the 1950s [7|. The cancellation of mass singular terms such as \n{m^/m e ) in the total 
rate, but not in distributions such as the electron energy spectrum, led to the development 
of the Kinoshita-Lee-Nauenberg theorem, which explains how to construct "infrared-safe" 
observables in quantum field theory where such effects cancel 0. The calculation of the full 
one-loop corrections in the SU(2) xU(l) theory of the electroweak interactions was one of the 
first such computations performed 0. The full two-loop corrections to the muon lifetime in 
the Fermi model, needed for a precision determination of Gf, were completed several years 
ago ; recently, the two-loop results in the full electroweak theory were obtained ll| . 

Muon decay continues to be of interest in particle physics. The TWIST experiment at 
TRIUMF l3| measures the electron energy and ang ular distributions in polarized muon 
decay; the first results were recently reported in [l3j. It is anticipated that TWIST will 
eventually measure the Michel parameters 14]|, which describe muon decay for the most 



general form of the four-fermion interaction, to a precision of ~ 10~ 4 . This significantly 
increases the sensitivity of muon decay to deviations arising from new physics. For exam- 
ple, the lower bound on the mass of the Wr in the manifest left-right symmetric model 
is improved from Mw R > 400 GeV to My/ R > 800 GeV, competitive with limits coming 
from the Tevatron, while the bounds on the left-right mixing parameter £ are improved 
by nearly an order of magnitude |5J). Such precision requires a careful consideration of the 
higher order corrections. As noted above, the radiative corrections to quantities such as the 
electron energy distribution contain large logarithms of the form ln(m^/m e ), which enhances 
their effect. The presence of mass singularities makes it impossible to compute the radiative 
corrections to the electron energy spectrum by neglecting the mass of the electron from 
the very beginning, the approximation that has been used successfully in the calculation of 
QED corrections to the muon lifetime [10]. This feature makes the calculation of the 0(a 2 ) 
corrections to the spectrum a challenging problem that has defied solution for many years. 

It was realized recently that the logarithmically enhanced parts of the second order QED 
corrections can be computed using the factorization of mass singularities traditionally dis- 
cussed in the context of QCD. In this way, the two-loop corrections with a double loga- 
rithmic enhancement, 0(a 2 ln 2 (m^/m e )), were calculated in (l5j |. and the singly-enhanced 
C(a 2 ln(m At /m e )) terms were computed in [16]. At the midpoint of the electron energy spec- 
trum, the sizes of these two terms are respectively —7 x 10~ 4 and 3 x 10 -4 . There are two 
interesting features of these results. The first is that the corrections are larger than the 
anticipated experimental precision, 10 -4 . The second is that the single-logarithmic terms 
are not a full factor of ln(m^/m e ) « 5 times smaller than the double-logarithmic terms, 
indicating that the naive power-counting based on the size of this logarithm might not hold. 
Both of these facts render a full calculation of the 0(a 2 ) corrections desirable. 
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In this paper we compute the 0(a 2 ) QED corrections to the electron energy spectrum 
in muon decay. The full dependence on the electron mass is retained. We use a method of 
performing next-to-next-to-leading-order (NNLO) calculations developed by us in a recent 



series of papers [17[. Our technique features an automated extraction and numerical can- 
cellation of divergences which appear as poles in the dimensional regularization parameter 
e = (4 — d)/2. In muon decay, ultraviolet divergences and divergences arising from soft 
photon emissions appear as 1 /e poles, while emission of photons along the electron direction 
is regulated by the finite mass of the electron and leads to logarithms of the ratio of the 
muon to electron mass, \a{m^/m e ). From the technical point of view, the fact that the 
electron mass plays the role of the collinear regulator leads to some differences as compared 
to calculations with only massless particles. Having masses as regulators reduces the com- 
plexity of the analytic structures which must be integrated over multi-particle phase-spaces 
or over virtual loop momenta. However, issues of numerical stability appear since multi- 
dimensional integrals are regulated by (m e /m^) 2 ~ 10~ 5 . We find that the presence of a 
mass regulator significantly simplifies the treatment of real emission processes of the type 
— > z/ M z/ e e + 77; however, the computation of virtual corrections becomes more complex, 
compared to a purely massless case. We describe these and other technical issues in detail 
in the main body of the paper. 

Many other physics applications require computations of higher order corrections to the 
decay spectra of massive particles. For example, the structure of the 0(a 2 ) corrections to 
muon decay is identical to the 0(a 2 ) QCD corrections to semi-leptonic b — > u and b — > c 
transitions, which are used to extract the CKM matrix elements \V u b\ and \ V c b\, the 6-quark 
mass and other important parameters of Heavy Quark Effective Theory The calculation 
of QED radiative corrections to the electron energy spectrum discussed in this paper can 
be easily extended to obtain differential results for semi-leptonic 6-decays at NNLO. In 
fact, some of the technical issues become simpler for 6-decays, particularly b —>■ c. In this 
case, collinear singularities are regulated by the factor (m c /mb) 2 ~ 4 x 10~ 2 , rather than 
(rrie/m^) 2 « 10 -5 , leading to more stable numerics. Precise predictions for heavy particle 
decay spectra will also be important for measurements at both the LHC and a future linear 
collider. Both experiments will search for anomalous top quark couplings through final- 
state distributions in its decay t — > bW, and will determine the CP properties of any scalar 
boson discovered through angular properties of such decays as ZZ, WW, ff |l9j . The 
techniques required to analyze higher-order corrections to these decay modes are very similar 
to those presented here. 

This paper is organized as follows. In the next Section we introduce our notation and 
discuss general aspects of the computation of the electron energy spectrum. In Section IIHI 
we describe our computation of the NNLO QED corrections. In Section HVl we discuss our 
results. Finally, we present our conclusions. 



II. NOTATION AND SETUP 

We discuss in this Section some basic notation needed to describe muon decay. We begin 
with the Lagrangian 

£ — £qed + £-f- (1) 
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£qed contains the kinetic terms for the fermions and photons, along with the QED inter- 
actions, 

Cqed = -\F, V F^ + ^ f [iP- m,\ ij; f , (2) 

/ 



while Cf contains the Fermi interaction, 



C F = -2V2G, 



(3) 



Here, P F 



75) /2 is the usual left-handed projection operator. The Fermi Lagrangian 



can be Fierz rearranged into the following form: 



(4) 



The QED corrections to this Lagrangian are finite to all orders in a |20f] after the fermion 
mass renormalization is included. Since the QED corrections do not affect the neutrino part 
of this Lagrangian, and experiments do not probe properties of the neutrinos, they can be 
integrated out to produce an effective ji — e current. We demonstrate this here. We first 
note that since QED interactions only affect the leftmost fermion bilinear of Eq. we can 
write the squared matrix element for the process /i — > ev e v^ + X as 



\M\ 



I KA pa I 



x Tr 



(5) 



Mf^e+x denotes the matrix element formed from the leftmost bilinear of Eq. 0] together 
with any QED corrections. Eq.(j5J) must be integrated over the appropriate phase-space to 
obtain the electron energy spectrum: 

dT 



— = j [dYl^ evv+x ]\M[ 



d P 2 Dt / [du^ epnt+x }\M p ;^ e+x \ 2 x / [du } 



|Tr 



i>v,P{p iv^laPh 



(6) 



where x = 2E/m fJL and E is the electron energy. In the second line, we have partitioned the 
phase-space so that the muon first decays into an electron, additional radiation denoted by 
X, and a massive state with momentum 



Put = Pp~Pe ~PX- 



(7) 



This massive state then decays into the muon and electron neutrinos. The neutrino portion 
of the phase-space can be integrated over to obtain 



[dU 



Tr 



with 



T. 



nt 



nplt 



Pnt,pPnt 



nt 

pa ! 



giving the expression for the decay spectrum 



9pa 



Pit 



dT 

dx 



dpl t / [dU 



}\M 



1 2n-\nt 
p— >e+X I pa' 



pa 



(8) 
(9) 

(10) 
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FIG. 1: A sample of LO and NLO diagrams which appear for the effective (j, — e current after 
the neutrinos are integrated out. The factor to be associated with the effective \i — e vertex after 
squaring the matrix elements is given in Eq. El 



Since, up to an overall numerical factor, T°* coincides with the polarization density matrix 
of a vector boson with mass p^ t , we can interpret Eq. El as the emission of such a boson in 
the fi — e transition. A sample of the diagrams that appear in this description are shown in 
Fig.HJ 

Since we are regulating divergences in dimensional regular izat ion, we briefly discuss our 
treatment of the 75 which appears in Pl in Eq. 0] Our conclusion is that a naive anti- 
commuting 75 can be used, and fermion traces containing an odd number of 75 matrices 
do not contribute. This follows from the observation that the tensor Tf^ is symmetric 
under p <-> cr; a contribution containing an odd number of 75 matrices will produce a form 
factor containing the completely anti-symmetric Levi-Cevita tensor and will vanish when 
contracted with . 

We now discuss the form in which we will present our results. We parameterize the 
electron energy spectrum with the variable x = 2E/m^ which lies in the range 



2m P a m 



2 



<x<l + ^. (11) 



We write the differential decay rate as 

^ = W £)>>(*). (12) 

dx 192vr 3 ^ \nj J K J K ' 

The LO and NLO results were computed in 0]; they can be obtained in a convenient form 
in 1(], 21 1 . The logarithmically enhanced contributions to f^ 2 \x) can be obtained from 0, 
161. The calculation of f^(x) beyond the logarithmic approximation, and including the full 
dependence on the electron mass, is the main subject of this paper. 



III. NNLO CORRECTIONS 

In this Section we discuss our computation of the NNLO corrections to the electron 
energy spectrum. We give a brief overview of the technical aspects of the calculation and 
then describe in detail the computation of double real emission, one-loop virtual corrections 
to a single photon emission, and two-loop virtual corrections. 
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A. Overview of NNLO corrections 



We first present a brief overview of the various components of the NNLO corrections. 
The differential decay rate contains a sum over several distinct components, 

where each dTy/dx is separately divergent and must be combined with the other components 
to produce a finite result. Our method of calculation follows the technique outlined in (l7| . 
We regulate both infrared and ultraviolet divergences in dimensional regularization, setting 
the space-time dimension d = 4 — 2e, and produce an expansion 

dx i=2 e l 

where the Aj are functions non-singular everywhere in phase-space. Since the A[ are non- 
singular, they can be computed numerically in four dimensions. The expressions for the 
dTy/dx can be combined, and the poles in e can be cancelled numerically. We must produce 
such an expansion for the following components. 

1. The real radiation corrections involve decays with two additional particles radiated 
into the final state. The two relevant processes are /i — > evv + 77 and fi — ► evv + 
e + e~. The first one begins at 1/e 2 , with the singularities coming from the phase-space 
regions where the photon energies vanish, while the second is finite. To handle these 
corrections, we use the techniques presented in [17]. We describe our phase-space 
parameterizations and discuss the extraction of singularities in Subsection 1111 Bl 

2. The real- virtual component includes the 1-loop virtual corrections to the process 
H — > evv + 7, and contributes beginning at 1/e 2 . We use a hybrid approach com- 
bining both analytical and numerical techniques to compute this component. We first 
use integration-by-parts identities and recurrence relations ji^] to reduce the correc- 
tions to a small set of master integrals. The recurrence relations are solved using the 
algorithm described in [23| and implemented in 24]|. We then solve for the master 



integrals numerically by applying the techniques of to their Feynman parame- 
ter representation. We discuss the details of this method, including how we handle 
imaginary components of the loop integrals, in Subsection IHI CI 

The virtual-virtual corrections contain the interference of two-loop virtual corrections 
to fi — > evv with tree-level diagrams. These begin at 1/e 2 . We deal with these 
completely numerically by applying the techniques of directly to their Feynman 
parameter representation. This numerical method of computing virtual corrections 
was pioneered in [2^, |26|. We apply it here for the first time in a fully realistic 
calculation, which includes tensor integrals and several mass scales. We discuss the 
details of this calculation in Subsection 1111 Dl 

We must include the square of the NLO virtual corrections, which contribute at 0(a 2 ) 
and produce poles beginning at 1/e 2 . Since the computation of this component can 
be performed with standard techniques, we do not discuss it further. 
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5. We must include both fermion mass renormalization, and external wave- function renor- 
malization. We renormalize in the on-shell scheme. The renormalization is performed 
by multiplying the LO and NLO results by the factor Zf x and by inserting the 
muon and electron mass counterterms into the NLO diagrams. Therefore, the mass 
counter-term is needed through O(o) only. For a fermion of mass m, the renormaliza- 



tion constants are 



Z 2 = 1 + £ 



n=l 



a T(l + e)m 
7r (47r)- £ 



-2c 



n 

7 {n) 
^2 1 



= _ 2 L_ Z (2) = 9 I 51 I 433 3 a( 3 ^ + _ l 3 ^ 

' 2 4e l-2e' 2 32e 2 64e 128 2 SV ; V ; 16 ' 



6m = ,n hm - m = rn(Z m - 1) = HSi + ^l^. (15) 

We note that contributions which arise from a closed fermion loop inserted into a 1-loop 
self-energy diagram, which appear at 0(a 2 ), have been removed from these formulae; 
they are more naturally included in the contribution discussed in the following item. 

6. Finally, we must include vacuum polarization corrections, in which a muon or elec- 
tron loop is inserted into a 1-loop diagram. These include the insertion of a closed 
fermion loop into a 1-loop vertex diagram, and insertions into external leg self-energy 
corrections, which contribute to the muon and electron wave function renormalization 
constants. These corrections form a finite subset. They can be computed using disper- 
sive techniques, as discussed in j3ll |32| , where their contribution to the muon lifetime 
and electron energy spectrum are computed. Since these corrections are discussed in 
the literature, and can be computed with standard techniques, we do not discuss them 
further. 

After combining items 1 — 6, both ultraviolet and infrared divergences cancel, leaving a finite 
result. This completes the brief summary of the terms which enter the 0(a 2 ) corrections; 
we now begin the technical discussion of their computation. 



B. Real radiation corrections 

We first discuss the contribution of the real radiation processes /1 — > evv + 77 and 
fi — > evv + e + e _ . A sample of the diagrams that contribute to these processes is shown in 
Fig. |21 They produce a contribution to the differential decay rate of the form 

= / d P l t J[dU^ ePnt+x ]\M^ e+x \ 2 T^. (16) 

In order to perform the integration in Eq. EE we must discuss both our phase-space pa- 
rameterizations and the singularity structure of the matrix elements. Denoting the radiation 
momenta by p 7 i, p 72 , p e -i, p e +2> we find the following denominators for each process: 

• n -»• evv + 77: = (p„ - p 7 i) - m 2 , c^ 2 , <P el = (p e + p yl ) - m 2 eJ dj 2 , c^ 12 = 
(Pa* - P7I _ P72) - m l, d 2i2 = (Pe + P 7 i + P72) - m l\ 
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FIG. 2: Sample diagrams which contribute to fj, — * evv + 77 (left) and fj, — > evv + e + e (right) 



— > ez/z/+e + e : <i 



12 



■ 2/ , ",2 - (/'. • /' 

(Pe+Pe-l+Pe+2^ 



e+2) , ^12 = {P^-Pe-l-Pe+2f- 



■m\, 



We first discuss our phase-space representation for the photon radiation process, which takes 
the form 



d P 2 nt 



[dU^ epnt+x ] 
2 



(2vr; 



3d-4 



ds nt J d d p nt d d p e d d p~ /1 d d p~ (2 



X <*(Pnt - S nt)5(Pe ~ ^D^l ) 5 (P^) 5 ^ <JV ~ ~ Pyl ~ Pf2 ~ Pnt) , (17) 



where the restriction of the electron energy fraction x is understood in the rightmost equa- 
tion. It is convenient to evaluate this in the rest frame of the muon and to choose the z-axis 
along the electron direction. In this frame, the momenta are 



p M = (rrv, 0, 0, 0) , Pe = {E e , 0, 0, /3E e ) , 
p 7 i = {Ex, Exsx, 0, Eicx) , p l2 = (E 2 , E 2 s 2 c^, E^s^, E 2 c 2 ) , 



(18) 



where E e , Ex, and E 2 denote energies, s\, s 2 , c\, and c 2 respectively denote sines and cosines 
of polar angles, and s^, denote the sine and cosine of the azimuthal angle. Following [l7l ]. 
we map this to the unit hypercube, and obtain 



AL f 1 dXtdX 2 dX 3 dX A dX 5 /% 2 2+2e (l + 5 2 
Jo ' 



X 



\4-4e, 



1-Ax 



,3-4e 



[Px 



l-2e 



where 

5 
E e 

iV 7 

«12 

/3 



x [A 2 (l - X 2 )] l - 2e [A 3 A 4 (1 - A 3 )(l - A 4 )P [A 5 (l - A 5 )f 1/2 - e 



m e /m M , ci = 2A 3 - 1, c 2 = 2A 4 - 1, = 2A 5 - 1, 



(19) 



x/2, 



A 2 (l-A!)(l + 5 2 -x) 



2 7 (27r) 



3d-4 



0„ 



2«i 
2vr d / 2 



£?2 



Ac 1 (l-A 2 )(l-A 1 )(l + 5 2 -a;) 



s nt = Ai(l + 5 2 



1 - /3c 2 ) 



r(d/2) 

A 2 (l-A 1 )(l + 5 2 -x)(l-n 1 -n 2 ) 



2/?i 2 

x), K\ — 1 



.2; 



(l-/5ci) 



l-45 2 /x 2 , rli • n 2 = c x c 2 + SiS 2 C0. 



(20) 



We have removed the integration over the energy fraction x, and have set the scale = 1; 
it can be restored with dimensional analysis. The matrix element denominators become 



d 1 



-A 2 (l-A 1 )(l + 5 2 -x) 



a [i2 



-/t 1 (l-A 2 )(l-A 1 )(l + 5 2 -x) 



«12 
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, 7 _ x\ 2 (l-\ 1 )(l + 5 2 -x)(l-0c 1 ) 



2m 

C2 



xmjl - A 2 )(l - Ai)(l + <5 2 - x)(l - ffca) 

2^12 



^ = _ (1 _ Xi){1 + p_ x) L + *a,(i - fa) + Wi - A a )(i - fa) 



el2 



1- Ax)(l + 5 2 -x) + — + 



2«i 

A 2 ki(1 - A 2 ) 



2k 



12 



«12 



(21) 



We note that all the divergences are produced by the overall multiplicative factors of (1— Ai), 
A 2 , and (1 — A 2 ); the bracketed terms in dZ 12 and cF el2 are finite for values of x away from its 
boundaries. In the language of |17[ , all singularities are factorizable; when the denominators 
are combined with the phase-space in Eq. EH the form of Eq. El can be produced by 
expanding in plus distributions: 



A 



-l+e 



-5(A) + E 



ra=0 



dXf(X) 



m"(A) 
A 

m w (A) 
A 



e 

n\ 



o 



■ dA /(A)-/(0), n , (A) 
A 



(22) 



We must now discuss our phase-space representation for the process fi — > evv + e + e , 
which takes the form 



dp 2 nt J [dU 



(2tt 



J ds nt J d d p nt d d p e d d p e - 1 d d p e+2 5(p\ 



nt S nt/ 



x S(p 2 e -m\)b(p\- x -m 2 e )S(p 2 e+2 -m 2 e )5 (d \p^ - p e -p e - 1 - p e+2 -p nt ). (23) 

It is convenient to view this decay as occurring iteratively; first, the muon decays into an 
electron and a massive "particle" with momentum p nt + p e -i + Pe+2- This massive particle 
then decays into p e -\ and another massive particle with momentum p nt + p e +2, which finally 
decays into p e + 2 and the neutrino pair. This motivates the following decomposition of the 
phase-space: 

ds nt ds nt i2ds nt 2 I1I2I3, 



(2< 



3d-4 



(24) 



where 



h 
h 
h 



d d p e d d p ntl2 S(p 2 e - m 2 e ) 5(p 2 ntl2 - s ntl2 ) <5 (d) (jV -p e - Pntu), 
d d p e - 1 d d p nt2 S(p 2 e - 1 - m 2 e )5(pl t2 - s nt2 ) S (d) (p ntl2 - p e -i - p nt2 ), 



J d d p e+2 d d p nt 5(p 2 e+2 - m 2 e ) 5(pl 



■Snt 



Pnt2 - Pe+2 ~ Pnt) 



(25) 



We evaluate I\, J 2 , and J3 in the rest frame of the massive "particle" that defines each 
phase-space. Doing so yields the following expression: 



d\id\ 2 d\^d\id\^ 



1 A34 Wl + 8 2 -X- 25) _,l-2e 



y/ S nt i2S n t2 

x [A 3 A 4 (1 - A 3 )(l - A 4 )P [Ab(1 - A 5 )]- 1/2 - e 



(3 6 xE 2 E i 



(26) 
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where 

^d-l^d-2^d-3 p _ g n t!2 + $ 2 — S nt 2 ,_, _ S nt 2 + S 2 — S nt 
25 + 4 e(27r) 3d-4 > 2- ^= , 3" ^= , 

1 + 5 2 - X, s nt = Ai (VI + 5 2 -x - 2(5) , s nt2 = A 34 A 2 + L 3A , 

2A 3 -1, c 3 = 2A 4 -l, c* = 2A 5 -l. (27) 

C2 is the polar angle of p e -\ within the frame of I2, while C3 and are respectively the polar 
and azimuthal angles of p e +2 within J3. To obtain the invariant masses that appear in the 
matrix elements, we must Lorentz transform the vectors p^, p e , p e -i, and p e +2 between the 
frames defined by the three phase-spaces. We note that all of the denominators that appear 
in the matrix elements are regulated by 5, and therefore the process p — > evv + e + e~ is 
finite. It is sufficient to set e = in Eq. EH] and to perform the numerical integration of the 
corresponding matrix element in four dimensions. 



N ee = 

S n tl2 = 

A34 = 

C 2 = 



C. Virtual corrections to single photon emission 



/1 e \i 




FIG. 3: Sample diagrams which contribute to p — > evv + 7 

Here we discuss the one-loop virtual correction to the process p — > evv+j; some diagrams 
that must be considered are shown in Fig- El This computation is conveniently performed by 
a combination of analytical and numerical methods. First, we express the integrals over the 
loop momenta through master integrals using the reduction algorithm of [2^ implemented 
in [2^. To give the list of master integrals, we introduce the following notation: 



DPi(ai, a 2 , a 3 , a 4 ) = / 



d d k 



(2vr) d (k 2 + 2p^k) a ^((k + p^ - Pl ) 2 - m 2 ^{k 2 + 2p e k) a *k 2a ± ' 

(28) 

Dp ( ) = f 1 

2lai ' ° 2 ' G3 ' G4j J (2n) d (k 2 + 2p fJi k) CL1 ((k + p e + Pl ) 2 - m 2 )^ {k 2 + 2p e k) a ^k 2 ^ ' 

We find that all the Feynman integrals needed for our purposes are expressed through fifteen 
master integrals that include the four-point functions DPi(l, 1, 1, 1) and DP 2 (1, 1, 1, 1), sev- 
eral three-point functions such as DPi(0, 1, 1, 1), DPi(l, 1, 1, 0), DP 2 (1, 1, 0, 1) and a number 
of two-point functions and tadpoles. These Feynman integrals depend on the energy of the 
external photon, c<j 7 ; when uo^ — > 0, some of the master integrals develop infrared singular- 
ities. The extraction of singularities is performed following |l7(; we Feynman-parameterize 
the master integrals, insert these expressions into our phase-space parameterization, and 
disentangle singularities in both the Feynman parameters and w 7 . 
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We write the real-virtual component of the NNLO corrections as 
dT RV 



dx 



J dpi J [dU^ epat+1 ]\M p ^ e+1 \ 2 T£ . (29) 



It is quite easy to construct a phase-space parameterization convenient for the extraction 
of singularities; it is very similar to the double real emission case discussed in the previous 
Subsection. We find 

J dpi J [dU^ epn J = 2 "y 2 "^_ 3 / dA^Ai- 2 ^ - A 2 r^ v (*, Xi, A 2 ), (30) 







where 

T (EL*z(l-zW-* 
rv (l-£ e (l-/3cosfl)) 2 ~ 2 ^ 1 ' 

Bmax = (l + 5 2 )/2, E e = x/2, z = E e /E max ,pl t = 2E max (l-z)(l-X l ), and costf = -1 + 2A 2 . 
In terms of these variables, the scalar products of the four-momenta s a b = 2p a ■ pb read 

s Me = 2E e , s M7 = 2cj 7 = s e7 = 2-E e w 7 (l - (3 cos 9). (32) 

1 — E e [l — p cos 6) 

From Eqs. EDI and E21 we see that potential singularities associated with the soft photon 
emission u; 7 — ► are factorized both in the phase-space and in the scalar products s a b'i 
therefore, their extraction proceeds along the lines described in (l7l |. 

An additional complication related to the real-virtual corrections is that some of the 
master integrals develop imaginary parts that, when the integral is Feynman-parameterized, 
appear as singularities in the integration region. This feature is very inconvenient since it 
makes it impossible to numerically integrate even otherwise finite expressions. We explain 
how we deal with this problem by considering the master integral DP 2 (1, 1, 1, 1) of Eq. |23 

We introduce a Feynman parameterization for the integral DP 2 (1, 1, 1, 1), and write it as 

DP 2 (l,l,l,l) = ^±^/dA 3 II Ml-Ex,)^, (33) 



i=1..3 i=l 



where 4> = — s e7 A3X 2 + x\ + (m 2 + s ei \^)x2 + xix 2 (s Me + X^s^). Changing variables to 
x 2 = AiA 2 and X\ — Ai(l — A 2 ), we find 

DP,(1, 1, 1,1) = /dA3dA ld A 2 — , (34) 

where A = (1 — A 2 ) 2 + (mg + s e7 A3)A2-l-A 2 (l — A 2 )(s Ate + A3S^ 7 ). The denominator in Eq. can 
become zero in the integration region, which makes accurate numerical evaluation impossible 
even for non-exceptional values of the photon energy. To circumvent this problem, we rewrite 
Eq. H2]in the following way. First, we integrate over Ai, producing a hypergeometric function 

DP 2 (1, 1,1,1) = / dA 3 dA 2 A 2 (- Se7 A 3 A 2 )— ^|f 21 ( 2 + e, 1 - e; 2 - e, 

(35) 
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We now use an identity that allows us to rewrite F 2 i(a,6, c, z) through i*2i(..., 1/z). We 
obtain 



' :r (2 + <) },< Jy w . , , -i / a V 2 " 



DP 2 (M, 1, 1) = / dA 3 dA 2 A 2 (- Sc1 A 3 A 2 , ]1+2e{ 



s e7 A 3 A 2 \ . T(l -e)r(l + 2e) ( A 



-1+e 



xf H 2 + f ' 1 + 2e ' 2 + 2f - a J + wU ) (36) 

For the hypergeometric function that appears in Eq. EE] we introduce a standard integral 
representation and arrive at 

-iT(2 + t) } ^ f A 2 Af 



dp 3 (i, i,i,i) = w n^ 



(4tt)^ 7 a 1 i (A - s e7 A 2 AiA 3 ) 2+e 

r(l -e)T(l + 2e) (-3 e7 A 3 A 2 ) ' 2 ' 
T(2 + e) s^AsA 1 -^ 



+ ~rvo ,a 7 \ ai f • (37) 



It is easy to see that the denominator of the first term on the right-hand side of Eq. E3 does 
not vanish inside the integration region; the imaginary part appears only from the second 
term on the right-hand side of Eq. which is oc (— l)~ 2e . Hence, Eq. ET71 can be used for 
numerical integration after disentangling singularities in u 1 and the Feynman parameters. 



D. Two- loop virtual corrections 

We compute the NNLO two-loop diagrams numerically. A basic ingredient of this ap- 
proach is the method of sector decomposition. In the past, this technique has been applied 
to scalar loop-integrals. In this paper, we extend the approach to compute a full two-loop 
amplitude; to the best of our knowledge, this is done here for the first time.. The tensor 
integrals that emerge in the two-loop amplitude can be expressed in terms of scalar integrals 
using, for example, the procedure in |29|, |30(. In principle, the latter could be computed 
with a brute force application of the algorithm in ji^, |26|. However, we have found that 
it is more efficient to adopt a slight modification of that approach. Specifically, since the 
number of Feynman diagrams we deal with is not large, we derive a Feynman integral repre- 
sentation for each of the diagrams. Such representations are not unique; the ones we derive 
simplify the evaluation of tensor integrals. First, we introduce a Feynman parameterization 
for the propagators of one of the loop integrals. We then integrate out the corresponding 
loop-momentum, and insert the result into the second loop. We carry out the remaining 
loop integration with a new set of Feynman parameters, using the approach of |25l . 126^ . 



li e n e fj, e 




FIG. 4: A sample of two-loop diagrams which contribute to fx — > evv. 
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As an example, we derive the parameterization for tensor integrals in the cross-triangle 
topology, which is the second diagram in Fig. HJ. We consider the integral 

r d d h d d k 2 {h} m {k 2 } n 

J m d l 2 A x A 2 A z A^As 1 ' 

where 

A x = kl A 2 = k 2 2} A 3 = kf + 2k x ■ p M , A 4 = k\ + 2k 2 ■ p e , 

A 5 = (k 1 + k 2 ) 2 + 2(k 1 + k 2 )-Pn, A 6 = fa + k 2 f + 2 (fci + fc a ) • p e , (39) 

and p 2 = m 2 = l,£>g = m 2 . We denote a tensor of rank m in the numerator with {k} m = 
^.A»i-Mm_ We g rs t introduce Feynman parameters for the propagators in the fc 2 loop. We 
write 

1 = r(4) Z 1 ^AitU 2 cU 3 A 3 (l - . v 

A 2 A A A 5 A 6 1 J 7o [(fc 2 + g) 2 -A 3 (l-A 3 )C a ] 4 ' 

with 

g = X 3 h + 77, 

r] = X 3 [Xip^ + (1 - Ai)p e ] + (1 - A 3 )A 2 p e , 

C a = fcj + 2k\ ■ p — — — — -, 

^ 3 U — ^ 3 J 

p = [X lV ^ + (1 - Ai - A 2 )p e ] . (41) 
We then shift the momentum k 2 , 

k 2 = K-q; (42) 
the shift yields a sum of tensors in K with ranks i < n: 

{M„-£*W*- (43) 
It is now straightforward to integrate out the loop-momentum K, using 

Wn = ( _d^ tM A ^r (44) 

mC z/ 2( ^2 + A)Q I ^ 2 «r(a) n ' 1 j 

T n = for odd n, and T = 1,T 2 = g^ 2 , T 4 = g w g i*si* + g nxm g ^ + g^g^^c. In 
order to perform the fci integration we introduce a new set of Feynman parameters A4, A5 
and shift the momentum k\. The shift yields a new set of terms which, after the integration, 
become 

Xil _r (2 + 2. - i±i W /' ( f[ ,A t ) - ^"^iL- , (45, 



\fc=l 



with 

= A 5 r/ 2 + A 3 (l - A 3 ) bvA 4 (l - A 5 ) + A 5 pf (46) 
The tensor integral is now written as 

x = X] fa • • • ' A 5 ; p m , Pe) Xij, (47) 
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where the terms are polynomials in the Feynman parameters; they are produced from 
the shifts of the loop momenta. The above Feynman representation is ideal for integrating 
over Feynman parameters after the sector decomposition in [3, 2(| is applied. Explicit 



expressions for the rather lengthy polynomials in Eq. 07| are not required in order to 
write down a Laurent expansion in e; these functions are only used in the Fortran code 
where coefficients of the e expansion are evaluated. We emphasize that our parameterization 
is very convenient since it allows us to treat tensor integrals on the same footing as scalar 
integrals. 

A technical complication arises when we consider two-loop diagrams with a self-energy 
insertion on either the muon or electron line, as in Fig. Quadratic singularities of the 
form A _2_e are produced, where A denotes one of the Feynman parameters. We cannot use 
the expansion of Eq. E21 to extract these singularities as a Laurent expansion in e. This 
occurs because one of the propagators in such diagrams appears squared. We solve this 
problem with the following procedure. As before, we start by performing one of the loop 
integrations; here it is easy to integrate out the self-energy. Let p be the momentum entering 

fj, e 




FIG. 5: A two-loop diagram with a self-energy insertion. 

the self-energy loop and m the mass in the propagator with momentum p. The result of the 
integration is 

8 = ( - 1) ' T(1 + e) /' dAl fi>Ar'(l-AQ- : , (48) 

£ Jo ( p 2 _ m 2 _ Aim^ 

where / is the polynomial from tensor reduction. We then insert this result into the second 
loop integration and obtain the following structure in the denominator 

A = tt^ r~p^- ( 49 ) 

(p 2 - m?f (p 2 - m 2 - ^g) 

A direct Feynman parameterization of A and sector decomposition produces quadratic sin- 
gularities. We avoid this by writing 

(hH^Y e ! x-l-e 

A = A^i— - e / d\ 2 ^ j-- (50) 

(P - m 2 ) 2 Jo (p2 _ m 2) ( p2 _ m 2 _ ^A^i-) + 

The first term in Eq. EDI leads to a straightforward one- loop integration, since all propagators 
are raised to integer powers. In the second term, the offending propagator is not raised to 
a quadratic power anymore. We employ a parameterization of this term following a similar 
procedure as in the cross-triangle topology discussed above. 

To summarize, we derive representations for the two-loop diagrams in the process fj, — > evv 
which (i) are amenable to sector decomposition, (ii) treat tensor and scalar integrals on the 
same footing, and (iii) are free of quadratic singularities. We then produce an e-expansion 
of the diagrams using Eq. 1221 and finally we evaluate the coefficients of the expansion 
numerically. 
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IV. RESULTS 



In this Section we give numerical results for the 0(a 2 ) corrections to the electron energy 
spectrum in unpolarized muon decay. We present our results in the form of a relative 
correction, 5^ = (a/7i) 2 f^ 2 \x)/f^(x), where and are respectively the LO and 
NNLO coefficient functions defined in Eq. This form allows us to study the magnitude of 
the corrections with respect to the relative experimental precision. We employ the numerical 
values = 105.658357 MeV and m e = 0.510998902 MeV and the on-shell value of the QED 
coupling constant, a = 1/137.0359895. For numerical estimates we consider electron ene rgie s 
in the range 0.3 < x < 0.95, which matches the acceptance of the TWIST experiment |l2j . 
We have checked that integrating our result over x reproduces the correction to the total 
decay rate found in (lOj, within numerical integration errors. 

As mentioned in the Introduction, the decay spectrum contains logarithms of the form 
In (m^/m e ), indicating that the electron energy is not physically observable as m e — > 0. The 
NNLO coefficient function can be expanded as a series in this logarithm: 



f^\x) = hi\mjm e ) ff>(x) + Hmjm e ) fT>(x) + tf>(x). 



e(2) 



(2), 



7(2), 



(51) 



The f^\x) and f[ 2 \x) terms have been calculated previously in [l5|, [l6{. The new result 



of this paper is the term ffi (x) 

the electron energy spectrum was previously estimated as ~ 10 
necessary to match the precision expected in the TWIST experiment. 



The uncertainty associated with the impact of f 



(2)/ \ 

"x) on 

4 , and its computation is 



o 




(2) 



FIG. 6: The ratio of the constant NNLO coefficient relative to the tree result, <5 L 
(a/7r) 2 /Q 2 ^(3;)//( )(x), versus the electron energy fraction x. The y-axis has been scaled by 10 4 . 



i2\ 

The magnitude of / i x ) relative to the tree-level result as a function of the electron en- 
ergy fraction x is presented in Fig. El To derive , we calculate using our numerical pro- 
gram, and subtract from it the logarithmically enhanced terms given in 2> 111 • We see that 



(2) 

for a large range of electron energies, the absolute value of Jq (x) is bounded by 0.5 x 10 



a value somewhat smaller than the theoretical expectations 15, 3]. To estimate the re- 
maining theoretical uncertainty on the electron spectrum, we note that 0(a 3 In 3 (m^/rrie)) 



15 



corrections to the spectrum have been computed in |33[; for moderate values of x, the cor- 
rections are in the range of few x 10 -6 . The pattern of logarithmic corrections at 0(a 2 ) 
indicates that the 0(a 3 In 2 (m^/m e )) terms might have a similar size. The hadronic cor- 
rection to the electron energy spectrum considered in [j^l is even smaller. Similarly, finite 
W^-mass effects are known, and influence the electron energy spectrum at the level of ~ 10 -6 . 
We take, conservatively, 5 x 10~ 6 as an estimate of the remaining theoretical uncertainty for 
values of x away from kinematic boundaries. 

10 



5 



O 

* — i 

x 

X 



-5 



-10 

0.4 0.6 0.8 

X 

FIG. 7: The ratio of the In 2 , In, and constant NNLO coefficients, as well as the total result, relative 
to the tree result, 5% = ( a / 7T ) 2 fx( x )l versus the electron energy fraction x. The y-axis 
has been scaled by 10 4 . 

An interesting feature of the electron energy spectrum is that effects of radiative correc- 
tions are enhanced by large logarithms of the ratio of the muon mass over the electron mass. 
The computation of the logarithmically enhanced terms in the spectrum is a much simpler 
problem than the full calculation reported in this paper. Since large logarithms are routinely 
exploited in theoretical physics for a simplified description, it is interesting to gain some ex- 
perience on how well this approximation works in various calculations. To do so, we compare 
the double- and single- logarithmic enhanced corrections with the full second order QED cor- 
rection to the electron energy spectrum in Fig. [7| We see that the ratio of the single logarith- 
mic term over the constant term obeys the expectation \fi (x)/ fo (x)\ ~ ln^/me) ~ 5, 
while the ratio of the double logarithmic term over the single logarithmic term doesn't: 
\f2 2 \ x )/ fi 2 \ x )\ < 5. Moreover, we note that the double-logarithmic terms overestimate 
the full correction. Because of the cancellation between the doubly- and singly-logarithmic 
enhanced terms, the relative importance of fo 2 \x) increases. For example, at x — 0.5, the 
constant term fo (x) changes the second order correction by about 10%; this should be 
compared with the naive estimate 1/ In 2 (m M /m e ) ~ 4%. From this we conclude that the 
leading logarithmic corrections give a correct order-of-magnitude estimate; however, the full 
result can deviate from the leading logarithmic approximation by a factor of 2 — 3. 
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V. CONCLUSIONS 



In this paper, we have presented a calculation of the 0(a 2 ) QED corrections to the 
electron energy spectrum in muon decay. The NNLO QED corrections, relative to the tree 
level result, are in the range —5 to 8 x 10~ 4 , depending on the electron energy. This is larger 
than the 10~ 4 precision expected from the TWIST experiment at TRIUMF. The corrections 
contain logarithmically enhanced terms of the form ln(m M /m e ), which have been calculated 
previously in |15l Il6j|. The new result derived here is the constant term without logarithmic 
enhancement, which influences the spectrum at the level of ~ 0.5 x 10~ 4 . The inclusion 
of this correction reduces the theoretical uncertainty below the anticipated experimental 
precision. We have argued the the remaining uncertainty is at the level of 5 x 10~ 6 , which 
is negligible for any foreseeable experiment. 

Although only the electron energy spectrum is considered in this paper, the computational 
method introduced is flexible enough to permit a computation of any distribution in muon 
decay, with arbitrary restrictions on the kinematic variables of the electrons and photons. 
The calculation reported here can therefore be extended in several ways. For the TWIST 
experiment, there are two natural extensions: (1) to include polarization of the muon, which 
is present in the experimental setup; (2) to also constrain the lab-frame angle cos 8 of the 
electron, in order to match the fiducial region used by the TWIST experiment in their first 
analysis, 0.50 < cos 9 < 0.84 0. 

There are several applications of our method beyond muon decay, where precise calcu- 
lations of the decay spectra of massive particles are required. Higher order corrections to 
semileptonic and radiative b decays are needed for extraction of CKM matrix elements and 
fundamental parameters in heavy quark physics, and in searches for new physics. Decay 
distributions of top quarks, Higgs bosons, and new massive particles will be precisely mea- 
sured at the LHC and at a future linear collider, and will be used to elucidate the underlying 
theory describing what is discovered. We anticipate that the techniques developed here will 
be useful in performing these analyses. 
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